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We present the theory of dynamical spin-response for the Kitaev honeycomb model, obtaining exact results 
for the structure factor (SF) in gapped and gapless. Abelian and non-Abelian quantum spin-liquid (QSL) phases. 
We also describe the advances in methodology necessary to compute these results. The structure factor shows 
signatures of spin-fractionalization into emergent quasiparticles - Majorana fermions and fluxes of Z 2 gauge 
field. In addition to a broad continuum from spin-fractionalization, we find sharp ((5-function) features in the 
response. These arise in two distinct ways: from excited states containing only (static) fluxes and no (mobile) 
fermions; and from excited states in which fermions are bound to fluxes. The SF is markedly different in Abelian 
and non-Abelian QSLs, and bound fermion-fiux composites appear only in the non-Abelian phase. 
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I. INTRODUCTION 

A time-lag of several millennia between the discoveries 
of ferromagnetic and antiferromagnetic (Neel) order, despite 
their great microscopic similarity, underscores the importance 
of the availability of experimental probes matching the phe¬ 
nomena in question. The lack of a characteristic macroscopic 
observable for the Neel state has an analogy today in the lack 
of any local ground state signatures of topological states of 
matter, for which the most natural diagnostics - entanglement 
entropy or topological degeneracies - are not readily accessi¬ 
ble to present experimental technology. 

The identification of clear-cut experimental signatures is all 
the more urgent - after a frustratingly long search following 
the original proposal of quantum spin-liquid states [1], there is 
no longer any shortage of theoretical models exhibiting ‘topo¬ 
logical’ quantum spin liquid states [2-4] ; and in the meantime, 
several frustrated magnetic materials have been identified as 
promising candidates to host QSL physics [5]. 

Perhaps the most natural local diagnostics for spin-liquidity 
involve the concomitant and characteristic fractionalised exci¬ 
tations above the featureless, long-range entangled, topologi¬ 
cally degenerate ground states. Due to the mismatch between 
the quantum numbers of such fractionalised excitations on one 
hand, and the selection rules for standard scattering probes on 
the other, experiments do not usually couple to a single frac¬ 
tionalised quasiparticle, instead exciting multiple quasiparti¬ 
cles, thereby producing a featureless continuum response. 

The dynamical spin response, which can be probed using 
conventional experimental techniques such as inelastic neu¬ 
tron scattering (INS) and electron spin resonance (ESR), is 
in principle sensitive to not only the ground state but also to 
a wide range of excited states, even at zero temperature. As 
such, it may in particular be sensitive to the presence of frac¬ 
tionalised excitations. Study of the dynamical spin response 
has proven to be fruitful in applications to one-dimensional 
systems, where it allowed one to establish a quantitative 
correspondence between the theoretically predicted correla¬ 
tion functions (obtained exactly using Bethe-ansatz), and the 


results of inelastic neutron scattering measurements [8, 9], 
thereby confirming the presence of S = 1/2 spinon excita¬ 
tions. 

In the recent Letter [10] we have commenced an analogous 
programme for a two-dimensional quantum spin liquid. The 
conclusions of Ref. [10, 11] were based on an exact calcu¬ 
lation of the dynamic structure factor for the celebrated 2D 
Kitaev honeycomb model (KHM); such exact results had thus 
far mainly been restricted to one dimension [12-15]. 

The Hamiltonian of the Kitaev model is remarkably sim¬ 
ple, having only nearest-neighbour exchange. This simplicity 
has led to a number of theoretical proposals for its realiza¬ 
tion in condensed matter, and in cold atomic systems [42, 43]. 
Materials whose spin and orbital degrees of freedom are 
strongly entangled in presence of spin-orbit couplings, such 
as {Na,Li} 2 lr 03 iridates [42, 44-48], and more recently a- 
RuCls [49-53], are currently the most promising candidates 
to realise Kitaev physics. Some of these are believed to be 
in the proximity of a quantum spin liquid state. Remarkably, 
residual high energy features of these putative QSLs might 
have already been observed in present systems [52-54], de¬ 
spite the fact that the latter are known to form a long-range 
ordered phase. 

The KHM represents one of the exceedingly rare instances 
of a tractable strongly-interacting quantum system in two spa¬ 
tial dimensions [16]. As a representative of a broader class 
of QSLs whose emergent degrees of freedom are Majorana 
fermions and Z 2 gauge fluxes, it has become an archetype for 
a QSL. Despite being formulated a decade ago, it still holds 
surprises, and is being actively studied, e.g. in the contexts 
of the calculations of ground state degeneracy [17], entangle¬ 
ment entropy [18], transitions between different topological 
phases [19-21], disorder effects [7, 22], global quench dy¬ 
namics [23,24], and the effects of doping [25-28]. There exist 
a number of integrable generalizations of the model [29-31], 
as well as its three-dimensional extensions [32-39]. 

While the calculation of the time-independent correlators 
is simple when expressed in appropriate variables [40], the 
calculation of dynamical correlators has turned out to be con- 
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siderably less straightforward. As noted already in Ref. [40], 
it is possible to map this calculation onto a non-equilibrium 
problem involving a quantum quench of a local potential, the 
physics which closely resembles the venerable X-ray edge 
singularity problem [41]. 

Here, we aim to provide a complete theory of the dynami¬ 
cal spin-response in two-dimensional Kitaev QSLs. We con¬ 
sider various different spin liquids - both gapped and gapless 
Abelian as well as gapped non-Abelian. The latter can ap¬ 
pear upon breaking time-reversal symmetry, and is of special 
interest due to proposals of using its non-Abelian excitations 
for topological quantum computations. For all of these, we 
provide the numerically exact dynamical structure factor, ex¬ 
tending our recent work in Ref. [10] 

We find a rich phenomenology, in which each of the consid¬ 
ered QSLs appears with distinctive signatures of the emergent 
fluxes and the Majorana fermions. Some of these properties 
are rather surprising, such as the appearance of a gap in the re¬ 
sponse for a gapless QSL, and the existence of a sharp (delta- 
function) response even for a fractionalised (both gapless and 
gapped) spin liquid. The explanation of the various features 
of the response are natural and simple in terms of the frac¬ 
tionalised degrees of freedom, e.g. involving the gap to a state 
with a pair of fluxes; or a bound state of the Majorana modes 
expected for a p-wave superconductor Hamiltonian represent¬ 
ing the non-Abelian QSL. It would seem hard even to ratio¬ 
nalise these phenomena in an alternative language. Therefore, 
this ensemble of results can be seen as a rather direct valida¬ 
tion of the fractionalised picture; given its richness, we do not 
provide a detailed list of our results in the introduction, and 
instead devote Section III to a non-technical account of our 
central results, which has been written with a reader in mind 
who is interested in phenomena but not too concerned about 
technical details. 

Finding the exact solutions presented here has led us to en¬ 
gage in a fair amount of method development. Much of the 
more technical material included here aims to give a reason¬ 
ably self-contained account of this. We have in fact devel¬ 
oped a number of complementary approaches, both exact (for 
finite systems based on determinant representation of correla¬ 
tion functions, and in the thermodynamic limit using singular 
integral equations) and approximate but simple (which we call 
the adiabatic approximation). It is perhaps worth noting that 
these should have applicability well beyond the present con¬ 
text. Much of this technical material has been collected into a 
set of appendices. 

While this paper presents an exact treatment of a particu¬ 
lar model QSL, its features should for the usual reasons be 
relevant to a much wider range of QSLs, qualitatively and 
(semi-)quantitatively. Namely, the gap in the response, which 
originates from the flux gap, the broad continuum due to spin- 
fractionalization, and the sharp delta-function response due 
to dynamical rearrangement of Majorana density of states, 
might hold for QSLs whose low-energy degrees of freedom 
are heavy fluxes of gauge-field coupled to dispersive fraction¬ 
alized excitations. These points are discussed as part of our 
closing outlook section. 

The structure of the remainder of this paper is the follow¬ 




Figure 1: (a) Honeycomb lattice showing bond directions x,y,z. 
Three-spin interactions (with coupling constant K) generate next- 
nearest-neighbour hopping for Majorana fermions along di. (b) Ma¬ 
jorana fermion dispersion at the isotropic point Jx = Jy = Jz, 
which is gapless for A = 0 (left) and gapped for A / 0 (right), (c) 
Phase diagram of the extended KHM. For A / 0 the ground state 
is a gapped non-Abelian QSL in the central triangle (grey), and a 
gapped Abelian QSL in the outer triangles (white). 

ing. In Section II we introduce 2D honeycomb Kitaev model, 
and its non-Abelian extension. We summarize our main find¬ 
ings for the dynamical spin-correlators together with the dis¬ 
cussion of their qualitative features in Section III. A brief out¬ 
line of Kitaev’s exact solution, is provided for completeness in 
Section IV. In Section V we present details of our calculations 
of the dynamic structure factor (SF). In Section VB we out¬ 
line two complementary exact methods for calculations of the 
dynamical correlation functions, and in Section V C a num¬ 
ber of approximate approaches. In Sections VIA and VIB we 
discuss qualitative features of the structure factor in the whole 
phase diagram of the extended Kitaev model. The main part 
of the paper follows with an outlook. Section VII, placing our 
work in a broader context and outlining directions for further 
work. 


II. EXTENDED KITAEV MODEL 

The KHM has spin-1/2 degrees of freedom on the sites of 
a honeycomb lattice. The spins interact via bond-dependent 
anisotropic Ising exchange J^, where the three directions la¬ 
beled by a = x^y^z distinguish the three bonds that share a 
given lattice site, as illustrated in Fig. 1. In the following we 
will also discuss an extended KHM, which is obtained from 
the original model by adding three-spin interactions. The lat¬ 
ter is generated by leading order terms in the perturbative ex¬ 
pansion in the strength of a small external magnetic field [16]. 
The three-spin interactions break time-reversal symmetry and 
generate a gap in the spectrum of Majorana fermions, giving 
rise to non-Abelian excitations [55]. 

The Hamiltonian of the extended KHM can be written in 
terms of Pauli matrices and we use the symbol {ij)a to in¬ 
dicate that two nearest neighbour sites (nn) i, j share the same 
a-bond. The extended KHM is obtained by adding next near¬ 
est neighbour (nnn) interactions between three spins 
associated with each pair of bonds {ij)a^{jk)b sharing the site 
j, where the direction of the component c is complementary 
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to a, 6. The Hamiltonian of the extended KHM model is 

^ = ( 1 ) 

nn nnn 

Ground states of the Hamiltonian Eq. (1) fall into three 
classes [16]. For iT = 0 there are two distinct phases, which 
are gapless and gapped Abelian quantum spin liquids. At non¬ 
zero K all of these phases acquire a gap, and the excitations of 
the formerly gapless state become non-Abelian. In all phases 
the independent degrees of freedom are static Z 2 gauge fluxes 
living on the plaquettes of the lattice, and dynamical Majorana 
fermions deflned on the sites. The time-evolution of the Ma- 
joranas is generated by the Hamiltonian whose form is flxed 
by a particular conflguration of the Z 2 gauge held. 

Our central objective is to calculate the dynamical structure 
factor (SF) 
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which can be measured in INS and ESR experiments. It is the 
Fourier transform in time and space of the dynamical correla¬ 
tion function 


St^{t) = (3) 

where af {t) is the a-th component of spin-operator in Heisen¬ 
berg representation at time t on site i. 

Here we consider zero temperature case, so that the aver¬ 
age (...) is taken in the ground state of the Kitaev model. The 
dynamical structure factor contains extensive information on 
excitations in the model, as is evident from the Fehmann rep¬ 
resentation 

- ^ 0 ]), (4) 

A 

where Ex is the energy corresponding to an eigenstate | A). 


energy difference between the ground state and the lowest- 
energy state with a flux pair in adjacent plaquettes. 

One component of the above-gap response is broad in en¬ 
ergy and results from Majorana fermion excitations. Its low- 
energy onset is at the two-flux gap in the phase with gapless 
Majorana excitations, but is higher in energy in the gapped 
Majorana phase. As we discuss below, this broad response 
is due mainly to either single-fermion or two-fermion excita¬ 
tions, depending on spin component and Hamiltonian param¬ 
eters. 

A further striking aspect of the dynamical response is that 
it includes in some instances sharp ((5-function) components 
in frequency space, a remarkable feature in view of the fact 
that all independent quasiparticle excitations of the model are 
fractionalised and so cannot be created individually by the ac¬ 
tion of a local operator such as af. In our discussion below, 
we identify two distinct physical mechanisms by which these 
sharp contributions arise. One mechanism involves ‘zero- 
fermion’ transitions, in which only Z 2 fluxes and no Majorana 
fermions are excited; the other stems from the bound states 
that are characteristic of vortices in this kind of non-Abelian 
spin liquid. 

While the excitation spectrum of the KHM is independent 
of the sign of exchange interactions, the ferromagnetic and 
antiferromagnetic models are clearly distinguished by the q- 
dependence of their response. Viewed in direct space, sign 
reversal for Ja leaves the on-site correlator unchanged but 
reverses the sign of the nearest neighbour In reciprocal 
space, this sign reversal transfers intensity in a characteristic 
way between the centre and boundary of the Brillouin zone, 
as examined in Sec. VIA. 

Much of the behaviour summarised in Fig. 2 can be under¬ 
stood starting from a selection rule for Majorana excitations. 
We set this out in Sec. Ill A, and provide a more detailed dis¬ 
cussion of our results for each phase in Sec. IIIB. 


A. Selection rules and a dynamical transition 


III. SUMMARY OF RESULTS 

We And that different representatives of the family of Ki¬ 
taev Hamiltonians encompass a set of qualitatively different 
responses, which we depict schematically in Fig. 2. All have 
in common that the spin correlations are ultra-short ranged, as 
first noted in Ref. [40] : the structure factor contains contribu¬ 
tions only from on-site and nearest-neighbour correlators, and 
only those with the same spin-component. This is a conse¬ 
quence of the static nature of the emergent Z 2 gauge fluxes. 
As a result, there are no sharp features in reciprocal space - 
in itself of course a classic ‘necessary but not sufficient’ diag¬ 
nostic of spin liquid behaviour - so that we restrict our plots 
in Fig. 2 to 5'q=o(cc). 

The next and considerably more surprising result is that, 
in all cases, the dynamical response is gapped, regardless of 
whether or not the underlying spin liquid phase has an exci¬ 
tation gap. The minimal gap in the response is given by the 


The states |A) that contribute to Lehmann expression, 
Eq. (4), for the dynamical structure factor are ones with non¬ 
zero matrix elements (A|(jj |0). Expressed in terms of Z 2 
fluxes and Majorana fermions, they obey selection rules which 
we now discuss. 

The flux selection rule is very simple (see also [40]): the 
action of cr^ inserts fluxes through the plaquettes either side 
of the bond {jk)a, as illustrated in Fig. 3. Since the ground 
state |0) is flux-free and fluxes are static, |A) contains this flux 
pair and no others. 

To introduce the selection rule involving Majorana 
fermions, consider in the first instance the ferromagnetic 
Abelian KHM deep in the gapped phase, with Jz ^ Jx, Jy > 
0 and AT = 0. It is then natural to discuss energy eigenstates 
in the basis of eigenstates of (jf, and for states in this basis 
to count the number of neighboring pairs of spins {ij)z 
that are antiparallel. One can define a spin parity operator 
Pz = (—1)^^ = Ylj The Kitaev Hamiltonian commutes 
with Pz and so all energy eigenstates can be chosen to have 
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(a) gapless Abelian (isotropic) 




(b) gapless Abelian (anisotropic) 



(e) gapped non-Abelian (anisotropic) 



(c) gapped Abelian 



Figure 2: (i) Schematic dependence of the structure factor 5'q=o(^) on energy uj and spin component a for the Abelian and non-Abelian 
QSL phases of the (extended) KHM. The panels (a) - (e) show behaviour at different representative points in the phase diagram along the line 
Jx=Jy, as indicated in (ii), with if = 0 for (a) - (c) and K > 0 for (d) and (e). Distinct spin components are denoted by solid (a = z) and 
dashed (a = x, y) lines. The insets show the density of states N{uj) of Majorana fermions and the energy of a Majorana bound state, where 
this is induced by presence of a flux pair. As discussed in Sec. Ill, sharp ((5-function) contributions to S^=q (ct;) appear for the Abelian model in 
the region of the phase diagram that is unshaded in (ii), but not in the shaded region. An additional sharp component is present in the extended 
KHM, throughout the non-Abelian phase and in some regions of the Abelian phases. 



Figure 3: (Colour online). A measurement of a dynamic structure 
factor leads to a sudden insertion of a pair of Z 2 gauge-fluxes (shown 
with minus signs in red). 


a definite parity. Moreover, is unchanged by the action of 
a single af operator, but is reversed by the action of a single 
(jf or operator. Hence the ground state |0) couples only 
to states I A) with the opposite parity for the structure factor 
components having a = x^y, and only to states | A) with same 
parity for the component with a = z. 

Deep in the gapped phase the lowest energy wave function 
for Majorana fermions in any flux sector has predominantly 
= 0, and therefore belongs to the = +1 sector. The 
higher energy states with single Majorana fermion excitations 
consist mostly of one antiparallel spin pair and belong to the 
Pz = — 1 sector. They form an energy band that is centred on 
2Jz and has a width set by Jx and Jy. 

Components of with a = x^y therefore arise in this 

part of the phase diagram only from states | A) that contain odd 
numbers of fermion excitations. Their main weight is due to 
single fermion excitations and is concentrated in a band near 


00 = 2Jz. These components also have some weight in higher 
bands near odd multiples of 2 J;^, but this turns out to be very 
small. 

Conversely, the component with a = z involves only states 
with an even number of matter fermion excitations. The low¬ 
est in energy of these is simply the unique excited state |Ao) 
containing no matter fermions and only the added flux pair. 
Since the matrix element (Aq |(3-J |0) involving this state is non¬ 
zero, it contributes to Sf^{uo) a (5-function with finite weight 
at frequency uo = Pao ~ ^o- Higher energy contributions 
form bands around even multiples of 2 J;^. 

We find a dynamical transition at which this ^-function in 
the structure factor disappears on moving through the phase 
diagram. The mechanism is as follows. Away from the limit 
Jz ^ Jx^Jy, the parities of the ground states in the relevant 
fiux sectors (zero-flux and the three two-fiux states that have 
fiux pairs either side of x, y or z bonds) are no longer nec¬ 
essarily even. In fact, their relative parities change on a line 
shown in Fig. 2(ii). All four ground states have the same par¬ 
ity for Jz ^ Jx^Jy^ but near = Jy the ground state 

with a fiux pair across a z-bond has opposite Pz parity to the 
other three ground states. As a result, there is no (5-function 
contribution to any component of Sf^{uo) in this phase. Note 
that the boundary on which this dynamical transition occurs is 
distinct from the previously-known thermodynamic boundary 
between the gapless and gapped phases, in fact lying within 
the gapless phase. 

Of course, similar arguments can be constructed using par¬ 
ity operators Px and Py based on the other components of 
spin, leading to the conclusion for the Abelian model that 
S^^{uo) has a ^-function contribution in a region of the phase 
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diagram where Ja is dominant, but not elsewhere. 

There is also a second mechanism that may generate a sharp 
contribution in the response. It is operative if the spectrum 
of Majorana fermion excitations that contributes to {|A)} in¬ 
cludes an isolated level, separate from continuous bands. Let 
|Ai) be the state containing only a flux pair and this fermion 
excitation. Provided (Ai|(j^|0) is flnite, this state contributes 
to a (^-function at frequency uj = Ex-^ — Eq. This is 

the case at some values of Jz ^ Jx = Jy in the extended 
KHM, for a = x^y within both dynamical phases, and for 
a = z in the dynamical phase that includes the isotropic point. 
Here the isolated level is a Majorana bound state trapped on 
the flux pair that is introduced into |0) by the action of It 
is known that flux excitations in the non-Abelian phase carry 
bound states of Majorana fermions, and these spatially local¬ 
ized states appear below the gap of the single-particle Majo¬ 
rana fermion continuum [ 22 ] . 


B. Qualitative features of the response 

We now discuss more fully the behaviour shown in Fig. 2, 
where we set JxjJz = Jy!Jz — 3 so that there are only two 
distinct components of the response: and = 

syy{uj). 

1. Abelian QSL 

Schematic results for the structure factor in Kitaev Abelian 
QSL phases are shown in panels (a) - (c) of Fig. 2 (i). 

In the isotropic model (a) is non-zero above the 

energy cost A for introduction of a flux pair. Its dominant 
weight arises from single Majorana fermion excitations, but 
a tail continues to higher energy. Although the energy width 
of the Majorana fermion band determines the extent of the 
main response, the energy-dependence of the intensity has no 
simple relation to the magnitude of the Majorana density of 
states, because the response involves fermion propagation in 
the presence of a flux pair. It nevertheless reflects features 
such as the van-Hove singularity. 

In the anisotropic model (b) there are distinct responses 
and including different flux gaps, and 

A;^. Within the gapless phase, both components have non¬ 
zero contributions above the respective flux gap. In addi¬ 
tion, beyond the dynamical phase boundary has a 5- 

function contribution at cc = A;^. 

In the gapped Abelian phase (c) the (5-function in at 

UJ = Az persists, but there is an energy gap separating it from 
the two-fermion continuum around uj = AJz. By contrast, 
S^^iuj) has no (5-function component and is dominated by a 
single-fermion continuum around uj = 2Jz. 


2. Non-Abelian QSL 

The response in the non-Abelian phase (AT 7 ^ 0 and 0.5 < 
j < 1 ) has features that are distinct from the ones which we 


And for the model with iT = 0. They arise because there are 
bound states of Majorana fermions associated with flux pairs. 

At the isotropic point (d) this composite flux-fermion bound 
state manifests itself as a single sharp component in the dy¬ 
namic structure factor, that would be absent from the corre¬ 
sponding Abelian phase. With anisotropy (e) the energies of 
sharp contributions to and are the sum of the fermion 
bound state energy and the two-flux gap, A^^ or A;^, and there¬ 
fore unequal. 

C. Broader implications 

The main features of the response described above are ro¬ 
bust against, for example, the addition of weak Heisenberg 
interactions to the Kitaev Hamiltonian, since spin-parity re¬ 
mains a good quantum number. The most important conse¬ 
quence of such additional interactions is that fluxes acquire 
dynamics. This will broaden the response around uj = A, but 
we expect that it will remain always gapped, and that distinct 
contributions to components of from states with zero, one 
and two matter fermion excitations will continue to be identi- 
flable. 


IV. REDUCTION OF THE SPIN HAMILTONIAN TO A 
FERMION QUADRATIC FORM 

The extended KHM model can be solved exactly following 
the original approach of Kitaev [16]. We introduce four Ma¬ 
jorana fermion species q and b^ with a = x^y^z on every 
lattice site i. These fermions obey the anti-commutation rela¬ 
tions {q, Cj} = 25ij and } = 2Sij5aa'- Spin operators 

can be represented in terms of Ci and b^ as 

= icih^. (5) 

Next, we deflne bond operators Uij = ib^b^ with i^j la¬ 
belling nearest neighbour sites at the ends of bond a. In 
terms of the bond operators Uij and the matter fermions c^, 
the Hamiltonian of the extended KHM reads 

H — ^ ^ iJdUijCiCj “h iK ^ ^ UijUjf^CiCj^. (6) 
{ij)a 

Bond operators are constants of motion with the eigenval¬ 
ues Uij = ±1. Thus the Hilbert space in which H acts can 
be decomposed into ‘gauge’ \E) and ‘matter’ \M) sectors. 
Replacing the bond operators by their eigenvalues we arrive 
at a Hamiltonian which is quadratic in Majorana fermions, 
and thus can be diagonalised. Note that three-spin interac¬ 
tions give rise to next-nearest-neighbour hopping for the mat¬ 
ter fermions. 

The Hamiltonian of Eq. ( 6 ) acts in an enlarged Hilbert 
space of four Majorana fermions at each site, rather then two- 
dimensional spin Hilbert space. This redundancy of the Ma¬ 
jorana mapping manifests itself in the local Z 2 gauge struc¬ 
ture, namely the physical properties (including the spectrum) 
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depend on the configurations {^o} of Z 2 fiuxes on the pla- 
quettes of the lattice, rather than configuration of bond vari¬ 
ables. The fiux on each hexagon is given by a product of 
bond variables = n(ij)eO physical eigenstates 

l^phys) = ^1^) are obtained using a projector to the physi¬ 
cal subspace P = ^P' [l + (—1)^^ (—1)^^]. Here P' is the 
sum of all operators which change bond fermion numbers in 
an inequivalent way [29, 56], and N^/f denote bond/matter 
fermion number operators. 

Observables should of course be evaluated using physical 
eigenstates |T^phys), but for the operators which do not change 
bond fermion number the same result can be obtained by 
omitting P and employing the unprojected states of the form 
|T^) = \F) (g) \M) (see [40] and Appendix A of Ref. [7]). In 
the following we restrict ourselves to observables of this type 
(note that for large systems complications from finite size ef¬ 
fects are negligible [57]). 

For a given configuration of bond variables {uij} the 
Hamiltonian can conveniently be written in the form 


yielding 

H = '^Enal^an - (13) 

n>0 n>0 

where En > 0 for n = 1... N are the eigenvalues, which 
depend on the fiux configuration, = F^n({0o})- The 
ground state of the matter fermion Hamiltonian Eq. (13) is 
defined by di\gs) = 0, with di = + yilfl- The 

ground state energy is therefore Egs = —^ 
der to find the global ground state of the spin Hamiltonian 
Eq. (1) one must compare ground state energies Egs{{(j)o}) 
in all fiux sectors. Fortunately, due to a theorem by Lieb, we 
know that the fermionic ground state in a translationally in¬ 
variant honeycomb lattice is fiux-free [59]. We denote the 
ground state of iT by |0) = |Eo) (g) |Mo), and fix the gauge 
suchthatii(^_^-)^|Fo) = +l|Fo) for all {ij)a. 

A. Ground state flux sector 


with the N X N matrix Mij = for N unit cells. Here 

caIcb is shorthand for the A^-component vectors CArIcBr- 
The next-nearest-neighbour matrices Fij and Dij vanish if 
AT = 0, but are non-zero at finite AT, see Eq. (6). We note 
that Eq. (7) is the most general form of a quadratic Majorana 
Hamiltonian. Instead of dealing with Majorana fermions it is 
more convenient to work with standard fermions, which can 
be obtained by combining two Majoranas into a single entity. 
To this end, we introduce two complex fermion species: bond 
fermions 

4).=® 

and matter fermions 

fr — '^{cAr (9) 

which obey standard anti-commutation relations [40]. Here 
A, B denote sublattice sites in the unit cell with coordinate r. 
The link variables Uij are simply related to the occupation 
numbers of bond fermions via 

Using the shorthand notation ca = P + f, CB=i{p-f) 
we write the Hamiltonian in terms of complex fermions as 



where h = (M M^) -h i{F — D) and A = (M^ — M) + 
i{F E D). The resulting Hamiltonian has the Bogoliubov-de 
Gennes form. It is diagonalized using a unitary transformation 
T, see e.g. Ref. [58], with TT^ = / and 

^ (a* A) 2” = (o P) <12) 


In the ground state fiux sector defined above, 
the Hamiltonian commutes with translations, and 
can be block-diagonalized via a Fourier transform 
~ SqGBZ ^ that 



In this representation Hq is equivalent to a BCS Hamiltonian 
describing a superconductor with a momentum-dependent gap 
Aq = — Hmsq — /^q (complex for 7 ^ 0 ), whose quasiparti¬ 
cle dispersion is ^q = ReSq, where Sq = 1 2 

and /^q = — 47f 3 5 sinqdi. Here ao = z^ai = 

X, 0(2 = y, the nearest neighbour vectors no = ( 0 , 0 ), ni = 
(1/2,73/2), na = (-1/2,73/2), and the six next-nearest 
neighbour vectors di,i = 1... 6 are defined in Fig. 1 (a). 

After writing the expression for the gap in the form Aq = 

I Aq I , the Hamiltonian Hq can be diagonalized by the Bo- 

goliubov transformation 



/ cos6>q sin 0q\ 

sin6>q COs6>q J 



(15) 


where 0^ is fixed by the condition tan 26>q = | Aq | /^q. Defin¬ 
ing Aq = ^q cos 26>q + | Aq| sin 26>q onc can write the Hamil¬ 
tonian in the form 


F0 = y],Bq(dldq-l/2), (16) 

q 


whose spectrum is given by 

Eq = 2^q+|Aq|2. (17) 

For AT = 0 the spectrum Aq = 215q | of fermionic matter 
excitations |Mo) is gapless if \Jz\ < \Jx\-\-\Jy\ (and per¬ 
mutations). At the isotropic point Jx = Jy = Jz there are two 
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Dirac cones positioned at Q = ±(27r/3, — 27r/3) with a linear 
energy spectrum ^(q) oc |q| at small energies, see Fig. 1 (b). 
In the presence of exchange anisotropy the Dirac cones move 
in the Brillouin zone, and merge at the transition line (between 
the gapped and gapless QSLs), so that for 17;^! > \Jx\ ^\Jy\ 
(and permutations) the spectrum is gapped. The phase dia¬ 
gram of the Kitaev model through the cut in the parameter 
space defined by -h -f = 1 is shown in Fig. 1 (c). 

The Dirac cones of the gapless phase (shown in grey) be¬ 
come gapped for nonzero K, see Fig. 1 (b). The spectrum 
remains gapless only along the dashed lines in Fig. 1 (c) with 
quadratic band touching at zero energy. The outer triangles of 
the phase diagram correspond to gapped Abelian QSLs whose 
fermionic bands are characterised by a zero Chern number. 
The particle/hole bands of the formerly gapless phase (central 
triangle) have Chern numbers = ±1, and the phase pos¬ 
sesses non-Abelian excitations [16]. 


V. CALCULATION OF THE DYNAMICAL STRUCTURE 
FACTOR 

In this Section we present several complementary methods 
which we developed to study the dynamical response in dif¬ 
ferent phases of the (extended) KHM. Those include two ex¬ 
act methods, as well as a number of approximate approaches. 
First, we outline the mapping to the quantum quench problem, 
which is the starting point of our analysis. Second, we present 
the exact determinant approach, which allows one to study nu¬ 
merically moderately large systems (see details in Appendix 
A), and may provide a starting point for further analytical in¬ 
vestigations. We also discuss another exact approach, based 
on the solution of an integral equation, which provides results 
in the thermodynamic limit. We conclude with the discussion 
of two approximations: the calculation of few-particle contri¬ 
butions to the response, and the adiabatic approach. 


A. Quantum quench correspondence 

The calculation of the dynamical spin-response in the (ex¬ 
tended) KHM can be mapped onto a local quantum quench 
problem for Majorana fermions, where a potential is created 
at time t = 0, and Majoranas propagate in the presence of 
this potential at f > 0, which is similar to a X-ray edge sin¬ 
gularity problem. This analogy was noticed by the authors of 
Ref. [40], and the results of our theory have been presented in 
earlier work Ref. [10]. Here we outline briefly the main steps 
of the quantum quench mapping; see Ref. [10] for details. 

A general expression for the dynamical structure factor is 
given by a Fourier transform of a two-time spin-correlation 
function. The latter can be expressed, following Kitaev, as the 
matrix element of Majorana fermions and fiux operators with 
respect to the ground state. In this representation the mapping 
to a quantum quench problem becomes very clear. The spin 
operator at time t = 0 in Eq. (2), for example on sublattice 
‘A given by af = ], contains bond fermion 


operators which exchange their number on the bond {ij)a be¬ 
tween zero and one. This corresponds to a change of signs of 
two fiuxes on two adjacent plaquettes sharing the bond, see 
Fig. 3. Since the fiuxes are static, the spin flip at time t in 
Eq. (2) must revert these fiuxes back to their original state in 
order to have a non-zero matrix element. This simple selec¬ 
tion rule leads to vanishing dynamical spin-correlators beyond 
nearest-neigbours [40], and the ones that remain non-zero are 
on-site, and nearest-neighbour correlators on the bond {ij)a 

S,f(t) (X (18) 

Elimination of the fiux degrees of freedom reduces the 
correlators to an essentially non-interacting form. How¬ 
ever, there is a price to pay, as in this representation one is 
faced with a non-equilibrium problem [40], whose physics is 
closely related to the celebrated X-ray edge singularity (see 
e.g. Ref. [60]). Explicitly, the ^-components of the correla¬ 
tors, which enter the structure factor read 

SABit) = 

SAAit) = (19) 

and similarly for the x, ^-components. The Hamiltonian Hz, 
which describes the time-evolution of Majorana fermions in 
the matter sector after the quench 

Hz = Ho + V, ( 20 ) 

differs from Hq only in the sign of the nearest- and (in the 
extended KHM) next-nearest neighbour Majorana hoppings, 
as can be seen from the form of a local ‘quench potential’ 
given by the sum of two contributions V = Vz ^Vk, where 

Vz = —‘^'iJzCACB, (21) 

Vk = 2iK[cA{cAd^ — CAde) + ^^(cBds ~ ^Bds)]- (22) 

For example the sign of the bond variable uab in Hz is oppo¬ 
site to the one in the ground state. 

This concludes the mapping of the problem of calculating 
dynamical spin-correlators in Kitaev model to a local poten¬ 
tial quantum quench. Despite an obvious similarity of the ex¬ 
pressions given above with the ones studied in the X-ray edge 
problem, we stress that the physics turns out to be quite dif¬ 
ferent, due to the presence of fractionalized quasiparticles. 

In the (extended) KHM the Majorana density of states at 
small energies either vanishes as zero energy is approached, 
due to Dirac dispersion, or has a gap, depending on the values 
of interaction constants Ja, AT. This low-energy behaviour is 
in contrast to what appears as an essential ingredient of an X- 
ray edge singularity problem, namely finite density of states, 
which lead to power-laws in the response. Related to this is 
the absence of the standard Anderson orthogonality catas¬ 
trophe [61] in our case. For example, deep in the gapless 
phase, the overlap between two Majorana ground states in dif¬ 
ferent fiux sectors (with and without V) does not vanish in the 
thermodynamic limit. There is another crucial ingredient in 
the Kitaev model, that is absent in the X-ray edge singularity 
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problem, and which leads to a new kind of Anderson orthog¬ 
onality catastrophe. Compared with the standard case, in the 
Kitaev model only the fermion parity, but not their number, is 
conserved. We find that this has a dramatic effect on the dy¬ 
namic correlation functions, and most remarkably, gives rise 
to a dynamical phase diagram: see Fig 7. 


B. Exact methods 

An exact evaluation of the dynamical structure factor start¬ 
ing from the Lehmann representation would amount to a sum¬ 
mation of infinite number of multi-particle processes gener¬ 
ated by a complete set of states |{Ac^}) = \Mq) in 

Eq. (34). Such a procedure is impractical, and instead we de¬ 
veloped two complementary exact approaches of calculating 
the SF, whose utility varies across the phase diagram. 


where is sl N x N diagonal matrix formed from the posi¬ 
tive eigenergies of the Hamiltonian Hz, and A’, 3^ are ma¬ 

trices correspond to a product of Bogoliubov transformations, 
see Appendix C, the matrix elements read 

= v/Det[A(i)][A-i(i)]gi. (27) 

Precise definition of the square root of the determinant can be 
found at the end of Appendix A. 

We note that one can use Eq. (27) in calculations for rela¬ 
tively large systems (we used a laptop to study systems with 
up to 10^ spins). In fact, in the gapped (extended) KHM 
phases it provides essentially numerically exact results for the 
response because finite size effects are negligibly small even 
in moderately sized systems. Remarkably, this method also 
works well near the isotropic point Jx = Jy = Jz because the 
time-dependent correlation function vanishes quickly with in¬ 
creasing time. This provided us with an independent check of 
the integral equation approach. 


1 . Determinant approach for correlation functions 


2. Integral equation approach 


Rewriting the correlators Eqs. (19) in terms of Bogoliubov 
quasiparticles dq which diagonalize the fiux-free Hamiltio- 
nian Hq [see e.g. Eq. (C2)] we obtain 

+ Yy)M{X*^ ± Fo*)]oo, (23) 

where plus/minus sign corresponds to AA/AB correlators. 
The main task is the evaluation of the matrix elements of the 
generic form 


Mgi{t) = (Mo|a,e-'"-*a;^|Mo), (24) 

where Hz is a Hermitian operator containing anomalous terms 
such SiS dqdk h.c. The latter conserve only the particle num¬ 
ber parity, but not their number. 

By representing the Eq. (24) in terms of a coherent state 
path integral one can obtain an expression for the matrix ele¬ 
ments Mqi in terms of Pfaffians 

Mgi{t) = ( 25 ) 

see details and definitions in Appendix A. 

This determinant approach is exact for finite-size systems 
in all phases of the (extended) KHM, namely it allows one 
to obtain the results for time-dependent correlation functions 
with any desired accuracy at arbitrary times. Similar ap¬ 
proaches are used in the studies of quantum quenches in 
e.g. quantum Ising model, non-equilibrium Luttinger liquids, 
as well as in the context of Full Counting Statistics (ECS) [62- 
64]. 

One can obtain a simplified expression for the matrix ele¬ 
ments which requires calculation of a single determinant, and 
an inverse of a A" x A" matrix at every time step of the calcu¬ 
lations, as shown in Appendix A. Here we only quote the final 
result. With the definition of a A x A matrix 

K = yle-^^^^y*F + xle^^^^XF, (26) 


In a recent Letter we showed that for a KHM it is possible 
to study the time-dependent correlators exactly in the thermo¬ 
dynamic limit [10]. Here, we outline the main steps of the 
calculations for completeness. 

In the interaction representation, with the time evolution 
governed by the Hamiltonian Aq, the local potential V plays 
the role of an interaction. The time evolution of an operator A 
in this representation has the form A(t) = ^ 


the wave-functions evolve under the ^'-matrix 


S{t,Q) = ^xex- 


;p[-z f c(tP(t)], (28) 

Jo 


where T is the usual time-ordering, and the nearest-neighbour 
dynamical correlator defined in Eq. (19) assumes the form 


SAsiO = -i{Mo\cA{t)S{t,0)cBmMo). (29) 


The main simplification and the reason why this mapping to 
an X-ray edge form of the correlator is possible can be traced 
back to a particularly simple local form of the impurity poten¬ 
tial e.g. for KHM Vz = —2iJzCACB, which is clear from the 
representation in terms of /-fermions. After introducing the 
occupation number operator fif = P f for the latter, where 
creates a complex matter fermion associated with the bond 
of the unit cell r = 0, the potential can be written as 


VPt) = -Uz[nf{t)-1/2]. (30) 

With these transformations the correlation functions can be 
reduced to simple expressions 


STB/AAit)=i{G{t,0)±G{0,t)], (31) 

where the two Greens functions are given in a standard time- 
ordered form 

G{t,0) = (32) 

G(0,i) = -i(T[/(0)/t(i)e-*/o‘drVAr)]), (33) 
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Figure 4: Cumulants (n = 0,1, 2) of Eq. (40) are shown for K = 
0 at the isotropic point. Black lines: exact result. Coloured lines: the 
single particle contribution. Note that, at the frequency of Majorana 
fermion band edge shifted by the flux gap (i.e. uj = QJz + A), the 
contribution from the single-particle excitations constitute 97.5% of 
the response. (Calculation is for a system with 65 x 65 unit cells. The 
small oscillations at low energies are due to finite size effects). 


These expressions for the GFs are similar to the ones which 
arise in the X-ray edge problem [41, 60], and can be evalu¬ 
ated exactly. This was done by us in Ref. [10] using a Dyson 
equation. The Dyson equation can be solved with the help 
of methods from the mathematical theory of singular integral 
equations, see e.g. a classic book by Muskhelishvili [65], and 
details of our calculations in the Supplementary Material of 
Ref. [10]. We checked that our numerical implementation 
of the determinant and integral equation approaches produce 
identical results (see also [66]). 

Note that deep inside the gapless phase, the ground states 
in the matter sector with and without a flux-pair have a flnite 
overlap in the thermodynamic limit, - remarkably, there is no 
Anderson orthogonality catastrophe [61]. This is due to the 
fact that the spectrum of matter fermions in the gapless phase 
is linear (Dirac-like), which leads to a vanishing DOS at small 
energies, and hence a small number of low-energy excitations 
which can be generated by an abrupt insertion of the fluxes 
(whereas in the standard X-ray edge problem the density of 
states is finite). Similarly absent are X-ray edge singularities 
in the response functions. 


3. Few-particle contributions to the response 

Before examining results from the exact solution, it is in¬ 
structive to look at the Lehmann representation of the dynam¬ 
ical correlation functions in the matter fermion sector. In the 
remainder we will concentrate on a discussion of the nearest- 
neighbour correlators S^^{t); the on-site correlators 
can be obtained in a similar way. First we deflne the basis | A) 
of many-body eigenstates of the Hamiltonian Hz with the cor¬ 
responding eigenvalues where Eq and Eq are the ground 

state energies of Hq and Hz respectively. After insertion of 


the identity operator \^){M i^^o Eq. (19) we obtain 
SAsit) = -i^e'*(^°-^x)(Mo|c^|A)(A|cB|Mo), (34) 

A 

whose Fourier transform gives 

= -2niY,{Mo\cA\X}{McB\Mo)d^_[E--Eo]- 

A 

(35) 

From this representation it is clear that the response van¬ 
ishes below A = Eq — Eq, which is the energy of the flux 
gap [10, 68, 69]. In a flxed gauge, Hq and Hz conserve mat¬ 
ter fermion parity and the only non-vanishing contributions to 
Eq. (35) arise from excited states |A) whose parity is oppo¬ 
site to the ground state \Mq). Therefore, the relative matter 
fermion parity of the ground states with and without fluxes 
plays a crucial role, and one has two possibilities. In case (I) 
the ground states of Hq and Hz have the same parity, in which 
case the states | A) must contain an odd number of excitations. 
In case (II), when the ground states have opposite parity, and 
I A) contains an even number of excitations. The sector with no 
excitations makes a special contribution in case (II), because it 
is just the ground state of the Hamiltonian Hz . Its contribution 
to is sharp in frequency, whereas the contributions 

from the sectors with flnite numbers of excitations are broad. 
As discussed in Sec. Ill, this striking distinction between cases 
I and II gives rise to a dynamical phase diagram because the 
relative parity varies as a function of coupling constants; see 
also Fig. 7. 

Using Eq. (34) one can obtain contributions from differ¬ 
ent fermion states e.g. ground state \Mf), single particle 
states |A) = h\\MF), two particle states |AA') = VJ)\,\Mf) 
etc., note the absence of tilde on A. 

In order to calculate these multiparticle contributions to the 
response one has to relate fermionic operators in the different 
flux sectors, as we explain in Appendix C. In case (II) (red 
region in the dynamical phase diagram) the approach must be 
adapted to the situation in which the relative parity is different; 
details can be found in Appendix D 1 , where it is shown that a 
sharp component, which appears exactly at the flux gap, arises 
from zero particle contribution in the Lehman expansion, 

cx5(w-A). (36) 

We And (see further discussion at the end of Sec. V C) that 
the single-particle contribution captures 97.5% of the total 
weight of the response at the isotropic point of the phase di¬ 
agram with FT = 0. Moreover, multi-particle contributions 
become smaller away from this point or at non-zero K, and 
so (depending on relative parities of zero-flux and two-flux 
ground states) single-fermion or zero- and two-fermion exci¬ 
tations account for nearly all the response. 












10 



AFM 



Figure 5: (Color online). The dynamic structure factor Ct;) for the ferromagnetic (left) and antiferromagnetic (right) KHM at the 

isotropic point, as a function of q and uj on the cut M—F—K—M through the Brillouin zone. 


C. Approximate methods 

1. Adiabatic approach 

Due to vanishing density of states in gapless KHM phases 
a replacement of an abrupt quench of the potential by an adia¬ 
batic switching on/off from — oo to +oo turns out to be a very 
good approximation [10] in the green region of the dynami¬ 
cal phase diagram (Fig. 7), and this replacement is exact in 
the low energy limit. In this adiabatic approach the potential 
generated by a flip of a flux-pair is switched slowly in time, 
thus the ^'-matrix 5(t, 0) can be replaced by 5(cxo, —oo), see 
e.g. [67], and the correlator assumes the form 

SAi\t) = -i(Mo|r[cA(t)cB(0)5(oo,-oo)]|Mo) 

=(37) 

with the only difference between Eqs. (37) and (19) being 
that the ground state |Mo) has been replaced by \Mf). This 
dramatically simplifles the calculations (because the integral 
equation can now be solved by a Fourier transform), and one 
arrives at an expression for the structure factor in this approx¬ 
imation 

( 38 ) 

A 

where X is a Bogoliubov matrix, see Appendix C. 


2. Lehmann representation and single-particle contributions 

Another quantitatively good approximation can be derived 
from the Lehmann representation Eq. (35). It holds in case (I) 
(green region in Fig. 7 (a)) where both ground states (with and 
without fluxes) have the same parity, so that only the states 
with odd numbers of fermion excitations parities contribute. 
Taking the |A) = b^lM^) we And the single-particle contri¬ 


bution to the structure factor 

= 87r|(Mf |Mo)|2 

(39) 

A 

see Appendix D 2 for details and deflnitions. Similarly, one 
can obtain two-particle contributions, see e.g. Fig. [12]. 

Comparison between cumulants 

nUJ 

dnn^Sq=o{n) ( 40 ) 

Jo 

of the exact solution with those of the single-particle response 
is presented in Fig. 4. The agreement is remarkable, which 
is due to the fact that only a small number of fermionic ex¬ 
citations is generated by the quench. The absence of or¬ 
thogonality catastrophe makes the single particle contribu¬ 
tion, Eq. (39), a quantitatively good approximation account¬ 
ing for 97.5% of the total intensity. Note that, above the 
single-particle Majorana band edge adjusted by the flux gap 
UJ = 6 J ;2 + A, the single-particle cumulants do not de¬ 
pend on cc. In contrast, the cumulants of the exact solution are 
frequency dependent due to many-particle processes (whose 
contribution to the response is very small). 

VI. RESULTS FOR THE STRUCTURE FACTOR 

In this Section we present a selection of the results for the 
structure factor in support of the schematic pictures shown in 
Fig. 2. The results here supplement those in Fig. 2 of our 
earlier paper [10]. 

A. Structure factor of the gapless Abelian QSL 

An overview of the wave-vector and frequency dependence 
of the structure factor, and a comparison between the ferro¬ 
magnetic and antiferromagnetic models, is shown in Fig. 5 
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Figure 6: Dependence of S^=o{uj) on u at the isotropic point of 
KHM. Blue line (highest peak) is a numerically exact result. Red line 
(middle peak), single particle contribution. Black line (lowest peak), 
the result of adiabatic approximation. The blue line is obtained in 
the thermodynamic limit via the integral equation approach; the red 
and the black lines are obtained for a system with 65 x 65 unit cells, 
with energy averaging over a window of width 0.025to remove 
finite size effects. Comparison between cumulants calculated within 
the same approaches is shown in the inset. 



Figure 8: S^=q{u) component of the dynamical structure factor in 
the gapless phase (green in Fig. 7) for Jz = 1. Red, green, and blues 
curves: f Jx — Jy — 0.8, 0.9,1.0. Note that the response diverges 
at the threshold on approach to the dynamical phase boundary. A 
similar divergence appears in the calculation of the the adiabatic ap- 
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Figure 7: (a) Dynamical phase diagram of the KHM model. The 
structure factor contains a (5-function contribution in the red 

region of panel (a), but not in the green region. Thermodynamic 
phase boundaries are indicated by white dashed lines, (b) Overlaps 
between Majorana fermion ground states from different fiux sectors 
on the line Jx — Jy (see Appendix. D 1 for notation). The weight of 
the (5-function, which is proportional to the overlap | |Mo) |^, is 
shown in the red region. At the dynamical phase transition, Jxj Jz — 
Jy!Jz ~ 0.71, the overlap drops to zero. The alternative overlap 
|(Mf|Mo)|^, shown in the green region, is finite where has 

no (5-function contribution. 



Figure 9: S^=o{uj) component of the dynamical structure factor in 
the gapless intermediate phase (red in Fig. 7), for J^ = Jy = 0.6, 
Jz = 1. Note a (5-function contribution to the response at the energy 
of the fiux-gap. The broad component of the structure factor shows 
significant multi-particle weight compared to other phases. 


for the isotropic point of the KHM. The main features of the 
uj dependence are relatively insensitive to q and to the sign 
of the exchange interactions. However, there is a striking 
shift of intensity from the center of the Brillouin zone in the 
ferromagnetic case to the edge of the Brillouin zone in the 
anti-ferromagnetic case, because of the change in sign of the 
nearest-neigbour correlator as discussed in Sec. III. 

The frequency dependence of cc) with q = 0 is 

shown in Fig. 6 for the ferromagnetic model at the same point 
in the phase diagram. The main features set out in Sec. Ill are 
apparent: response is non-zero only above the two-flux energy 
gap, and the dominant contribution extends over the energy 


width of the Majorana fermion band. Fig. 6 also demonstrates 
that single Majorana fermion excitations account for the ma¬ 
jority of the response, and that our adiabatic approximation 
[scaled using the sum rule Eq. (El)] captures the behaviour 
quite accurately. 

Evidence for the dynamical phase transition, which is dis¬ 
cussed in detail in Sec. Ill, is presented in Eig. 7. The structure 
factor includes a (5-function contribution in the indi¬ 

cated region of the phase diagram, with a weight that drops 
discontinuously to zero at the transition. On approaching the 
transition from the opposite side, a broad peak above the flux 
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Figure 10: Spectral features on the line j—JxlJz—JylJz for K — 
O.IJ^, shown in order to highlight the signatures of the thermody¬ 
namic phase transition between Abelian and non-Abelian phases at 
j = 0.5 and the dynamical phase transition at j ~ 0.73. Red line: 
overlap |(Mo|M^^)|; blue line: the flux gap A. Both the overlap 
and the flux gap are continuous at j = 0.5, but their derivatives are 
discontinuous, signalling the thermodynamic transition [21], while 
at j ~ 0.73 the overlap drops abruptly to zero, indicating the dy¬ 
namical phase transition. Black crosses: dependence of the energy 
of the lowest excited state (the band edge) for a translationally in¬ 
variant system (black crosses); green diamonds: energy of the lowest 
excited state for a system with two fluxes (green diamonds). For the 
values of j where these energies are different, the fluxes support a 
static (as opposed to a dynamically generated) fermion bound state. 


gap (see Fig. 6) sharpens and shifts to lower energies, see 
Fig. 8. The continuum response has a divergence at the flux 
gap energy precisely at the transition. After crossing the tran¬ 
sition, this sharp peak splits into a (5-function, and a flnite con¬ 
tinuum response as shown in Fig. 9. 


B. Structure factor of the extended KHM (K / 0) 

At a flnite value of the three-spin interaction K, the next- 
nearest neighbour hopping gaps out the Dirac cones (see 
Fig. 1), and the resulting QSL state hosts non-Abelian excita¬ 
tions. This can be seen in the representation of the Hamilton- 
ain in terms of spinless fermions, Eq. (14). These interactions 
break time-reversal symmetry, and the “superconducting gap” 
in Eq. (14) acquires a phase. The Hamiltonian in the transla¬ 
tionally invariant system (in the absence of fluxes) is identical 
to a Hamiltonian describing px + ipy superconductor, which 
can support bound states in vortex cores. More speciflcally, 
an isolated half-vortex, equivalent to a Z 2 flux in the KHM, 
carries a zero energy state [16, 70, 72]. 

When evaluating the dynamical structure factor we are con¬ 
cerned with fermion states in the presence of a pair of Z 2 
fluxes, induced by the action of aj on the flux-free found state. 
Since the two fluxes are in adjacent plaquettes, the zero energy 
modes that would be associated with each one if they were iso¬ 
lated are hybridised, forming a more conventional bound state 
at flnite energy, or merging with the continuum for some val- 
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Figure 11: Dynamic structure factor S^{uj) at Jx=Jy=Jz and 
K = O.lJz as calculated from the exact determinant approach 
(blue). In addition we show the single-particle contribution (red) 
and the adiabatic approximation (black, with total intensity flxed 
by the sum rule). Inset: Majorana fermion density of states N{uo) 
of the flux-free sector. Note a (5-function contribution to AS'q“(ct;) 
from the localized Majorana state bound to a flux pair. Its energy, 
A + Ef — 0.9247^, is below the single particle gap. The calcu¬ 
lations are done for a system of 56 x 56 unit cells, using an energy 
broadening of 0.025to reduce flnite-size effects. 


ues of Ja. The behaviour of this bound state and some other 
features of the Majorana fermion spectrum in the extended 
KHM are illustrated in Fig. 10. 


C. Fermionic bound state signatures in the structure factor 

In Fig. 11 we present results for the dynamic structure fac¬ 
tor at the isotropic point in the non-Abelian phase. The main 
features of the response are a (5-function component due to a 
Majorana fermion level bound to a flux pair, and a broad com¬ 
ponent from excitations to the fermion continuum. The sharp 
feature is at an energy cc = A + E[ = 0.9247^ which is the 
sum of the two-flux gap A = Eq — Eq — O.MbJz and the 
fermion bound state energy E[ , while the onset of the broad 
feature is at A + E 2 , where E 2 = 2.0027;^ is the energy of 
the fermion band edge. 

In Fig. 12 we show the dynamic structure factor at 

a point in the Abelian phase of the extended KHM. It displays 
some of the main features introduced in Sec. Ill: a distinc¬ 
tive difference between { 00 ), dominated by single-fermion 

excitations, and which has a sharp zero-fermion con¬ 

tribution at the energy A = 0.09057;^ of the flux-gap, and a 
small two-fermion band, with an onset at A + 2E[, where 
E[ is the lowest fermion excitation energy, lying at the band 
edge since there is no fermion bound state for this choice of 
interaction strengths. 

We examine in Fig. 13 the evolution across the phase dia¬ 
gram of the energy and weight of the (5-function contribution 
to in the extended KHM. This arises through different 

mechanisms in the two dynamical phases. For j < 0.73 it is 
due to a transition to a state with a flux pair but no fermion 
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Figure 12: Dynamic structure factor for inequivalent compo¬ 

nents X in the Abelian anisotropic QSL at Jx—Jy—^-2^Jz and 
K = O.IJ^. Inset: Majorana fermion density of states N{uj) in flux- 
free sector. The zz-correlator has a (5-function contribution at the 
two-flux energy, arising from a ‘zero-fermion’ excitation, and also 
a weak two-fermion contribution. The main contribution to the xx- 
correlator is from single fermion excitations. Calculations are done 
for a system of size 56 x 56 unit cells using a broadening 0.025 
to reduce flnite size effects. 



■Ix^fz fy/fz 


Figure 13: Dependence of the weight of (5-function contribution to 
on j — JxjJz — Jy/Jz^oxK = 0.1Jz. Inset: dependence 
of the energy Ct;min of the (5-function contribution and the two-flux 
gap on j. The thermodynamic and dynamical transitions occur at 
j = 0.5 and j = 0.73 respectively. Weights and energies vary con¬ 
tinuously across the transitions, although the origin of the (5-function 
is different in the two dynamical phases. 


excitations, and is visible only in the component a = z. For 
0.73 < j < 1 it arises from a transition to a state with a flux 
pair and a fermion excitation in a bound state, and appears in 
all components a. It is notable that evolution is continuous 
across the dynamical transition and also across the thermody¬ 
namic transition at j = 0.5. 


VII. OUTLOOK 

We have presented a complete and exact calculation of the 
dynamical structure factor of the various quantum spin liq¬ 
uid phases of the Kitaev honeycomb model, with results sum¬ 
marised in Section III. These contain both general features in¬ 
dicating the presence of fractionalised excitations, as well as a 
sufficient amount of detail reflecting the particular QSL phase 
to be useful as identifying diagnostics. 

A question which follows almost reflexively concerns the 
applicability of such a set of results to models away from spe¬ 
cial points of exactitude, or indeed actual materials. A case 
in point is the kind of Heisenberg-Kitaev Ji — J2 — Js — ^ 
model which has been adduced to account for properties of 
the Ir-Kitaev materials in 2D. 

Here, we expect important gross features to be robust, as 
the mechanisms underpinning these phenomena are not pred¬ 
icated on integrability. For instance, the emergent fluxes will 
continue to be natural variables to describe the spin liquid, 
even if they are no longer immobile or entirely absent from the 
ground state. A suppression of low-energy scattering due to 
the selection rule on the fluxes should therefore persist. Simi¬ 
larly, for the sharp delta-function features, both the ingredients 
underpinning their appearance survive departure from integra¬ 
bility. Given these two - parity selection rules, and the pres¬ 
ence of Majorana bound states in the case of the non-Abelian 
QSL - can further be distinguished by considering different 
components of the structure factor in the anisotropic spin re¬ 
sponse, this may be a particularly promising way for detecting 
a phase with non-Abelian anyons. Regarding the latter, how¬ 
ever, it is worth bearing in mind that there is the possibility 
of the bound state energy being increased to the extent that it 
will merge with the single-particle continuum. 

Our results indicate that the dynamical structure factor al¬ 
lows for a quite a detailed level of Majorana spectroscopy, 
reflecting e.g. their bandwidth, their interactions with the flux 
pair from zero- all the way to many-particle signals, their ar¬ 
rangement in the perturbed ground state or the presence of 
bound states. We believe that this is an important feature in 
the broader quest for Majorana physics, which has been cen¬ 
tral to topological condensed matter physics for a while. It 
will remain to be seen to what extent some of the flner fea¬ 
tures, such as the higher multi-particle continua with their rel¬ 
atively small weight, will in practice be visible. In the first 
instance, both the broad and sharp features at relatively low 
energies are going to be the most likely signatures, and our 
exact solution will hopefully be of use in modelling these in 
an attempt to fit experimental results [53]. 

The technology developed here to study the dynamical 
structure factor may be applied to other members of the large 
and growing family of Kitaev QSLs [29, 30, 36, 38, 39]. Ma¬ 
jorana spectroscopy as outlined above will likely be an ap¬ 
propriate framework for interpreting the results in most cases, 
where it is an entirely open question how ‘details’ - such 
as spatial dimensionality, dimensionality of the zero-energy 
Majorana manifold, symmetries and possibly different selec¬ 
tion rules - will manifest themselves. Beyond building up 
a compendium of possible behaviours, we hope that such a 
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programme will lead to a practically useful field guide for the 
identification of such QSLs. Considering the properties of sur¬ 
face states of 3D QSLs appears to be another promising line 
of future work. 

Finally, the general methodological advances of our work 
should be applicable in contexts beyond the KHM. An ap¬ 
proximate approach to the computation of the dynamic struc¬ 
ture factor, here called the adiabatic approximation, has been 
known and used elsewhere for a long time, even in the remark¬ 
ably similar problem of a missing core electron in a ’single 
graphite plane’ [78]. 

Beyond this, however, our calculation of the dynamic struc¬ 
ture factor can be taken as an exact solution of a local quan¬ 
tum quench [40], related but not identical, to the X-ray edge 
problem [41], and as such represents a contribution to non¬ 
equilibrium quantum dynamics in its own right, for which we 
have developed two complementary approaches. First, in or¬ 
der to obtain exact results in the thermodynamic limit we have 
adapted a method from the theory of singular integral equa¬ 
tions [10, 65, 71], which can be extended to calculations of 
the full frequency dependence for the local impurity quenches 
with non-standard fermionic density of states. Second, we 
have derived an exact determinant expression for an arbitrary 
quadratic quench (allowing for the presence of anomalous 
terms), which can be useful for numerical calculations with 
systems considerably larger than those accessible to exact di- 
agonalization. 

We hope that the insights presented in our work, together 
with the technology developed here, will be of use for the in¬ 
vestigation of an extended class of fermionic many-problems. 


VIII. ACKNOWLEDGEMENTS 


We thank G. Baskaran and F. H. L. Essler for helpful dis¬ 
cussions. The work of J.K. is supported by a Fellowship 
within the Postdoc-Program of the German Academic Ex¬ 
change Service (DAAD). J.T.C. is supported in part by EP- 
SRC Grant No. EP/I032487/1, D.K. is supported by EPSRC 
Grant No. EP/M007928/1. This collaboration was supported 
by the Helmholtz Virtual Institute “New States of Matter and 
their Excitations”. 


Appendix A: Determinant approach for correlation functions 


Below we present the derivation of exact expressions for 
the dynamical spin correlators, similar in spirit to those of full¬ 
counting statistics (FCS). These results can be applied to eval¬ 
uate numerically the structure factor in finite systems, whose 
size is much larger than is amenable to exact diagonalization 
approaches. The results of this Section are not constrained to 
a specific model, and may provide a starting point for further 
analytical considerations. 


1. Definitions 


The dynamical spin correlation functions of the (extended) 
KHM can be expressed in terms of the matrix elements, see 
Eq. (23), which have the following form 

Mgi{t) = (Mo|a,e-*"^*aJ|Mo), (Al) 


where Hz = Hq Vz is the Majorana (matter) Hamiltonian 
in the presence of two dipped fiuxes (playing a role of a local 
potential for Majoranas), and |Mo) is the ground state of Hq, 
defined as ag|Mo) = 0 for all q G BZ. 

We will now derive the expressions for the matrix elements 
(Al), which are suitable for numerical evaluation, in terms of 
Pfaffians. Eirst, it is convenient to express the Hamiltonian 
Hz in terms of operators d^. In this representation the Hamil¬ 
tonian assumes the Bogoliubov-de Gennes form 




hijalaj + + ^Aijaja] 


(A2) 


where d, A are the matrices to be defined later (the following 
discussion does not rely on a specific form of the matrices, 
provided that the Hamiltonian is Hermitian). 

The matrix elements defined in Eq. (Al) can be calculated 
using Grassmann path integrals. Eirst, we represent the state 
I Mo) in terms of a trace over complete set of states {|n)} using 
a projector to the ground state, which has no fermions, so that 


= y^(n|ajva]y_j ...a\ a^e di... aAr|n) 

n 

= Tr{e“*^^^djdi... dArd]^ ... djdg}. (A3) 

Here N is the total number of momentum states in the Bril- 
louin zone (which is equal to the number of unit cells). After 
commuting all aj to the right and all dq to the left, we obtain 

Mqi{t) = (—l)‘^^^Tr{djY ... dj^^dj_^... d| 

X ... dq-idq^i ... dAr} 

- SqiTi{d^^ ... dje-'^^'di... div}. (A4) 


Notice that a single creation operator dj , and a single annihi¬ 
lation operator dq is absent in the first term. Now our task is 
to derive a generating functional F[J*, J]. Then the matrix 
elements Mqi (and other Green functions) are obtained by a 
differentiation with respect to the source terms J. Eor exam¬ 
ple 


Tr{d]y...d|e ... dAr} 


d^^F[J*,J] 


Using the identity operator for the set of Grassmann vari¬ 
ables {0a}, a = 1.. .N (see e.g. Ref. [75]) 


N 

a=l 


d(f) cx^djcf) Q,^ 


-Ea 




(A5) 
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one can write for Vq = Tr[e 


-iHzt] 


Vo = 


f ^ 

/n 

^ a=l 




—iHzt I 


matrices of the Bogoliubov transformations have been defined 
in Sec. IV, and have a form 


f-i 




and T = 


A'* 3^* 
3^ A’ 


(All) 


Using Suzuki-Trotter decomposition on the r.h.s. of this equa¬ 
tion by inserting M — 1 times the identity operator, so that we 
have M time intervals = tjM), we obtain 

Po= lim 

M^oo J 


with A', 3^ being N x N matrices. We recall that this trans¬ 
formation relates the operators do, and ba, which diagonalize 
the Hamiltonian Hq in the fiux-free sector, and the Hamilto¬ 
nian Hz in the two-fiux sector respectively. Now we can write 
Eq. (A8) in a diagonal form. After introducing another vec¬ 
tor combined of Grassman variables (^^1 = = 

we obtain 


Here /£>[</)*, (^] = / [I^i na=i and the fields 

obey anti-periodic boundary conditions = —0^. Next, 
we introduce a Fourier transform in the Matsubara space 

P=f-i 

= E (A 6 ) 

p=-^ 

with the Matsubara frequencies given by ^ (p + ^) • 

The sum over frequencies in the Eq. (A6) can be separated 
into positive and negative components (note that co’_(p+i) = 
—ujp). We absorb the common factor i6t into a redefinition of 
matrices h = ih6t and A = iA6t, and take the large M limit. 
Note that the latter introduces a phase ambiguity in the path 
integral, see e.g. Ref. [73], which will be resolved at a later 
stage. By introducing a vector notation (in the index a), we 
combine the Grassmann variables into a single vector |T>^) = 
[ 1) 10“ ^^ ^, and after defining the matrix 


r -Er=o(^^l 

Vo= 

N N 

= JJ = JJ (1 + *). (A12) 

n=l n=l 

The last line was obtained by evaluating the Grassmann 
integrals, see e.g. [73], and we used a standard formula 
In {iujp + = ln[l -h (note that = i5tE^). 

2. Generating functional 

We construct a generating functional by adding source 
terms, represented by vectors of Grassman variables \J'p) = 
—I to the path integral 


lUJnj 


n 


-n 


i^p) 


j~ p — 


iujp E 

. ^At 


- 

I o ' ^ 


(A7) F[J*,J] = J 


the expression for the Vq assumes a compact form 




(A13) 


/ M-i 

In order to evaluate the path integral in Eq. (A8), we first di¬ 
agonalize the quadratic Hamiltonian 



h 

A ■ 

dj 


-'H. 

kJ 


= Ha 


(A9) 


The Gaussian integrals are calculated using Bogoliubov trans¬ 
formation, and after taking the inverse Fourier transform back 
from the Matsubara frequency space 


M 




jm 

^ r\j 


Vm 


/A ' 

m=l 


we arrive at the expression 


(A14) 


using a Bogoliubov transformation described by a matrix T, 
see Ref. [58], such that 


THf-^ 


n 0 

0 -Cl 


(AlO) 


where ClissiN x N diagonal matrix of eigenvalues Note 
that the eigenvalues in Eq. (12) and the eigenvalues 
in Eq. (AlO) differ by a factor of idt. We can now write the 
expression for H in the diagonal form H = — 

\ X]n>o In the remainder we will omit the constant con¬ 
tribution, which will be restored in the final expression. The 


E 


F[J*,J] =Poe 


~ 0-i^p{rn-n) 
iujp 

0 


0 


^+iujp{m — n) 
xuj p ^2 


f\jn) 


In the following we will only require the functional deriva¬ 
tives taken at the times having the index M, e.g. the matrix 
elements of the form d‘^¥/dJ^^dJ^*, because all operators 

d|, dj in Tr{... d|... ... dj ...} are taken at the same 

time (at the boundary, see Eq. (A6)), and we set m = n = M. 
Further we extend the sums over p negative frequencies, and 
introduce the definitions = n^{E), 
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where ^ is s vector of eigenvalues En. The phase ambigu¬ 
ity in the path integral can be fixed by comparing the results 
with the ones calculated by standard operator theory of low- 
order matrix elements, e.g. TT{dle~^^^^dq}, see Ref. [ 66 ]. 
The functions 


n^{E) 


1 

1 + 


(A15) 


assume the form of the Fermi-distribution, and we obtain 


where the sign factor (— 1 )^ on the second line appears due 
to anti-periodic boundary conditions on Grassmann variables. 
Differentiation with respect to Grassmann variables is essen¬ 
tially the same as integration (see Ref. [74]) so that 

L{t) = J dJN... dJidJ*...dJ*^ F[J*,J]. (A 20 ) 

After relabelling [Jl ■ ■ ■ J%Ji ■ ■ ■ Jn] —> [^i---^ 2 iv] we arrive 
at the result 


F[r,J]=Voe 


n-{E) 

0 


0 

n+(E) 


f\j) 


(A16) 


Note the factor of 1/2, which is due to the fact that the sums 
have been extended to negative frequencies; we have dropped 
the index M. The matrix can now be written in explicit form 


~n-{E) 0 

T — 

~A B' 

0 n+{E)_ 

± — 

C V 


where the entries S, C, D are given by x matrices 


Ar(Ar-l) f InTc-n 

L{t) = {-l)^—Vo X / d 02 N ■■■deie^‘^ 

^ ^ N(N-l) 

= (A21) 

In the last line we used Gaussian integration of anti-symmetric 
matrices with Grassmann variables, which results in a Pfaf- 
fian, see e.g. Ref. [74]. 

The first term in the matrix element Mqi in Eq. (A4) has the 
form with a creation and an annihilation operator removed, 
namely 


A = X'^n-X* +y'‘n+y, B = X^n-y* xy'fn+X, 
c = y'^n-X* +X^n+y, V = y'^n-y* + X^n+X. 


In terms of these matrices, the partition functions reads 

F[J*, J] = 

(A17) 

where we substituted V with —A^, which fixes the phase am¬ 
biguity. Note that in the absence of the anomalous terms in the 
Hamiltonian, the expression for the generating functional that 
we obtained reduces to the standard free-fermion result [75]. 

One could now in principle calculate the matrix elements 
in Eq. (A4) by differentiating Eq. (A 17) with respect to the 
sources. However we find it more convenient to proceed in a 
different way. First, we reorder the expression under the ex¬ 
ponent into a trace with the anti-symmetric matrix S = —S^, 
defined as 


-B A' 

-A^ -C ’ 


so that we can write 


(A18) 


(A19) 


Rgi (t) = Tr{a]v ... ...a[ 

X ... dq-idq+i ... ajv} 

where the measure of the integral is defined as 

D[0qi] = d02N • • • d02N-l+ld02N-l-l • • • dOq^idOq-i ... dOi. 

In the series expansion of the exponent, the only non¬ 
vanishing terms which remain after the integration are those 
which contain 2N — 2 Grassmann operators 

/ n — N-\-l r 

2— N-\-l ^ 

^ (VTTljj • • • = Pf'5[gZ], 

where P is a transposition of 2N — 2 indices (i.e. 1... 2N 
excluding q and 2N — 1. The matrix S^qij is obtained from the 
matrix S by removing two lines and two columns at positions 
2N — I and q. Finally we arrive at the exact expression for the 
matrix element 


3. Matrix Elements 


.V .X (Ar-i)(Ar-2) 

Rgiit) = (-1)-^- Vo (A22) 


Let us now discuss the calculation of the matrix element 
L{t) = Tr{d]y ... ... d^} 

= j D[4>*... 4>% 


The Eq. (A22) requires evaluation of a Pfaffian of the ma¬ 
trix with size 2{N — 1) x 2{N — 1). If numerical imple¬ 
mentation of this equation the determinant |Po| can become 
very large, while the absolute value of Pfaffians very small, 
which would bring large numerical errors. In order to regu¬ 
larise the matrix elements one can use the well-known prop¬ 
erty of Pfaffians, namely Ff[BAB^] = DetP x PfA for any 
matrix B of the same size as A. To this end we introduce a 
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2N X 2N diagonal matrix B = diag[(l + e )’ 

ATtimes 

whose determinant is equal to 

DetB = 11(1 + ) = T>o. (A23) 

n 

This matrix provides a required regularisation of the matrix 
elements through the expression 

VoFfS = Ff[BSB]. (A24) 

While the matrix under the Pfaffian on the l.h.s. of this expres¬ 
sion can become ill-conditioned, the r.h.s. of this expression 
can be evaluated numerically without difficulties. 

For the Pfaffian entering expression for the Rqi{t), we use 
the expansion formula of Eq. (B2) to derive the identity 


Pf5[,,] = (-l)^+^(-l) 


\A^+l 


Pf 


■ 0 0 ■ 

0 0 0 -fl 0 

0 0 

0-1000 
0 0 


, (A25) 


=5 


Ui} 


where the 2N x 2N matrix on the r.h.s. is obtained from 
the matrix S by setting the rows and columns 2N — I and q to 
zero, and Sq 2N-i = — 1 and S 2 N -1 g = +1- We have 

R^i{t) = (-l)«+*(-l)^^PoPf[5{5i}], (A26) 

where the Pfaffian is calculated for a 2A/' x 2N matrix. 

As was shown above for the matrix 5, one can regularise 
the Pfaffian FfS^qiy using the same matrix B. After defining 
a phase-factor 


K{t) = (A27) 

and collecting all contributions, we obtain the expression for 
the complete matrix element 

= K{t){Pf[BS^giyB] - 5giPi[BSB]}, (A28) 


where a precise definition of the square root is the following 

\/Det[A(t)] = v/|Det[A(i)]|e'‘"^(*)/2^ (A32) 

and the phase (pA{t) = arg[Det[A(t)]] is taken to be a con¬ 
tinuous function of time. Note also that Mqi{0) = 6qi. The 
phase contains a large linear part —2Eot, and in numer¬ 

ical calculations it is convenient to separate this contribution 
first. In fact the loop contribution in the diagrammatic expan¬ 
sion can be written as 

{S{t, 0)) = = e*®°VDet[A(t)]. (A33) 


Appendix B: Some useful properties of Pfaffians 


The results of the previous Section have been expressed 
in terms of Pfaffians, which often appear in Gaussian inte¬ 
grals over anti-commuting variables, see e.g. [74]. Below we 
present a short overview of the definitions and the properties 
of Pfaffians that are relevant for our calculations. 

A Pfaffian is an extension of a determinant for skew- 
symmetric matrices A = —A^. It is always possible to write 
a determinant of a skew-symmetric matrix as the square of a 
polynomial in the matrix elements, [74] 

[PfA]2 = DetA. (Bl) 


The formal definition of a Pfaffian for a 2N x 2N skew- 
symmetric matrix A is 


PfA 


1 

2^N\ 


®Sn[P] 

PEii...i2N 


^i2N-li2N 


with the matrix elements aij, and sgn[P] = ±1 the sign of 
the permutation P. Hence, the Pa ffian has the unique sign for 
the square root Pf A = ±A/DetA. Note that the Pfaffian of an 
odd-dimensional matrix is zero. 

Several properties known from determinants carry over in a 
modified way to Pfaffians [74, 76]: 

• Multiplication of a row and a column on a constant is 
the same as multiplication of the entire Pfaffian on this 
constant. 


which is Eq. (25) of the main text. This expression can be 
further simplified by generalising a theorem for Pfaffians, see 
Ref. [77], and we obtain 

Mqi {t) = -K{t)VoPfS X [Im + G\i , (A29) 


• Interchanging two rows and the corresponding columns 
flips the sign of the Pfaffian, 

• Adding multiples of a row and the corresponding col¬ 
umn to another row and a column does not affect a Pfaf¬ 
fian. 


where G is the upper right N x N block of the inverse of S, 
and /tv is the A" X identity matrix. One can still reduce the 
size of the matrices which needs to be inverted by a factor of 
two due to a special structure of the result. Let us introduce a 
N X N matrix 

A = (A30) 

where is a A x A diagonal matrix of . We arrive at 
the following simple result 

Mgiit) = VDet[A(i)][A-i(t)],,, (A31) 


Another very useful property (in particular for numerical com¬ 
putations) is an expansion formula for the Pfaffians 

2N 

PfA = ^(-l)*auPf[Aii], (B2) 

i=2 

with the reduced matrix An having the first row and the i-th 
column removed. In addition, we will use the relation 

Pf[BAB'^] = Pf A DetB, (B3) 

which is valid for an arbitrary square matrix B having the 
same dimension as A. 
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Appendix C: Ground state overlap from the Bogoliubov 
transformations 


The Bogoliubov transformations for the matter fermion op¬ 
erators in a given flux sector (in a flxed gauge) read [58] 

h =Xlak + Ylal (Cl) 

f] =Y]iai+Xla], 

where XX are x matrices, and we assumed summa¬ 
tion over repeating indices. The calculation of the dynami¬ 
cal structure factor requires matrix elements between states 
in different flux sectors. Let h and d be the operators in in¬ 
equivalent flux sectors in which Hamiltonians for the matter 
fermions is diagonal. For deflniteness let’s take a system with 
(index F) and without (index 0) fluxes. The Hamiltonian in 
each case is diagonalised by respective Bogoliubov transfor¬ 
mations, namely in the flux-free sector we have 


Xo* 

Yo Xo) 

with the inverse 

W 4 / 

In the sector with fluxes 








(C2) 


(C3) 


(C4) 


Now one can use these to express b operators in terms of d’s 


A’* 3^* 
3^ A’ 




(C5) 


Here we introduced the matrices 


A’* 

y X 


/ V* J_ V* vt _L V* 

[YFX^ + XpY^f YfY^+XfXIJ-^ 


The ground states are related by the unitary transformation 
T of Eq. (C6), and their relative parity is even or odd if the 
determinant of real orthogonal matrix 


B = UTU^ and C/ = . (C7) 

is equal to ±1. Provided that the ground states of Majoranas 
in the sectors with and without fluxes have the same parity, the 
two-flux ground state \Mf) deflned as h\MF) = 0, is related 
to the flux-free ground state d|Mo) = 0 via the following 
equation, see Ref. [58] 

\Mf) = Det (C8) 

Tij = (C9) 

Hence the overlap between two ground states reads 

{Mf\Mo) = Det(XtX)V (CIO) 

A related overlap arises in the X-ray edge problem between 

electron ground states with and without a local potential. 


1. Zero- and two-particle contributions 

In order to calculate the ^-function contribution in case (II) 
one has to modify slightly Eq. (34). The previous Appendix 
C explained how to relate different flux sectors via proper Bo¬ 
goliubov rotations, in particular Eq. (C8) establishes the rela¬ 
tion between ground states of different flux sectors. Impor¬ 
tantly, the product of the parities of two ground states can be 
found from Eq. (C7). It is clear from Eq. (C8) that \Mf) gen¬ 
erated using a proper Bogolyubov rotation has the same parity 
as I Mo), so that for case (II) the approach has to be extended 
to the situation in which the parity changes, which we do in 
the following. 

The problem with the naive use of the Lehmann representa¬ 
tion can be traced back to the fact that we have not projected 
the identity operator to the physical subspace. In order to cure 
this problem one can re-introduce projectors, or alternatively 
use improper Bogoliubov transformations. However, there is a 
simpler way, which is to take advantage of the gauge structure 
of the Kitaev model. As we discussed above, the model con¬ 
serves parity of the total number of fermions N = Nf + N^. 
A gauge transformation changes the parity of bond and matter 
fermions while keeping the total parity intact. We note that 
only the relative parity of matter fermion ground states in two 
flux sectors (which differ by local fluxes) is important. Since 
Majorana fermions are their own adjoints c\ = 1, a correct 
form of the Lehmann expansion can be obtained by plugging 
a modifled identity operator ca ^q. (19), 

which gives for case (II) 

STsit) = -ze'®°‘^(Mo|e-*"-*|A)(A|c^CB|Mo). (Dl) 

A 

Here we used the fact that CAC = e with 

Hxy = caHzCa = Hq F Vx F Vy (D2) 

is a gauge transformation, i.e. it does not alter the flux sector. 
The Hamiltonian Hxy can be obtained from Hz by revers¬ 
ing signs of all nearest neighbour, and next-nearest neighbou- 
rur hoppings on bonds sharing a site A. The eigenstates of 
Hxy form a set of many-body eigenstates |A), which can be 
generated by creating excitations on top of the ground state 
with fluxes, which it is convenient to represent in terms of 
the ground state without fluxes via Eq. (C8). Note that the 
spectrum is invariant under the gauge-transformation dis¬ 
cussed above, but the parity of its ground state \M^) is op¬ 
posite to the parity of the ground state of Hz namely \Mf). 
Hence, in the case (II), i.e. with vanishing overlap (M^lMo) 
due to different parities, the overlap of the gauge-transformed 
ground state is flnite, and is given by 

\{MF\Mo)\^ = 'Det\X^y\, (D3) 

where the matrix is deflned in the Appendix A. 
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The zero-particle contribution is obtained by restricting the 
sum over A in Eq. (Dl) to the ground state 

= - 2[Y^(Y - J'*X*)]oo}, 

SO that the structure factor at q = 0 reads 

= 87r|(M^^|Mo)|"<5(:^ - A) 

X {1 - [ytyjoo _ Re[rt^*x*]oo}. (D4) 

In order to satisfy the parity constraint, the next non-vanishing 
contribution to Eq. (Dl), in addition to the (^-function, arises 
from two-particle excited states 

= 87r|(M^^|Mo)|" 

xY,S{lo-[E^ + E^,+A]) 

XX' 

x{|GAA'|" + Re[iGtA'Gi*v]/2} (D5) 
with matrix elements 

4 % = {Mo\b\b\,\Mp^) 

= {Mpy\Mo) {yxix;^, + yxiJ^ikXlx' }, ( d 6 ) 



and 

41 = {Mo\cBoCAob\b{,\M^py) 

= “'*41 + flil + 9x1' ■ (D7) 

Here we also defined the following contributions 


Figure 14: Dependence of the nearest-neighbour equal-time correla¬ 
tor Sab {t — 0) on the values of the exchange couplings. Panel (a) 
corresponds to the case K = 0, where the strength of the n.n cor¬ 
relator is shown along the cut Jx + Jy + Jz = 1 in the parameter 
space. Panel (b), nearest-neighbour correlator as a function of K/ Jz, 
at several fixed values of j — Jxj Jz — Jyj Jz, see inset. 


5^1 = 2i{Ml^\My){Xx^i^Xx'jY^^ 

- XxiYioXx'jXjo + Y^YoiyxiXfy}, (D8) 

and 

5il = 2i{Mp<\Mo){yxiXE.,X^jEjiYio 

+ yxiXuX^XyjYjo + yxiXfiY,oXyjXjo+ 

^xiXioyyifFijYjo - YxiYioyx'iTijXjQ 

-yxiXiky^x'yo^jY;,}. (D9) 


2. Single particle contributions 

In case (I), green region in Eig. 7, the ground states of Ma- 
jorana fermions have the same parity in both (initial/final) fiux 
sectors. In order to calculate the spin-correlators from the 
Lehmann representation we insert the identity operator into 
Eq. (19) such that only the states of opposite parities con¬ 
tribute. Here we restrict the calculation to account for the 


contributions from single particle states. The latter have the 
form I A) = ). With the help of the equation 

Xxj-VxiEij=[X^]l^ (DIO) 

the nearest-neighbour spin correlator assumes the form 

A 

X {X -Y)lAix\'^~%{X + Y)jo, (Dll) 

where we use the summation convention on repeating indices. 
Erom this expression (and a similar one for 5'^^ ^) the single 
particle contribution at q = 0 follows, Eq. (39). 

Appendix E: Static correlators, and sum rules 

As a check of our calculations of the dynamical correlation 
functions we make use of the sum rules, e.g. 

s^^{t = 0) = ^ r duj (El) 

J — oo 
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which relates equal time correlators to the integrated dynami¬ 
cal response. We employ this as a check of our exact results, 
as well as the quality of multi-particle approximation. The 
equal time correlator can be obtained from the equation 

5yB(i = 0) = 2 ^ cos0q (E2) 

qGBZ 

with 6>q defined in Eq. (15). At the isotropic point for K = 
0, Ja = 1 we obtain = 0) = 0.5249 in the thermody¬ 


namic limit. In Fig. 14 (a) the equal time correlation function 
is shown in the full phase diagram [40] for if = 0. The iso¬ 
lated Ising dimer limit 5'^^ (t = 0) = 1 is quickly approached 
for Jx^Jy ^ Jz^Jz = 1- In panel (b) of the same figure we 
show the evolution of S^^{t = 0) as a function of K for 
different values of j = JxjJz — JyjJz^ Static correlations 
always decrease with increasing \K\. 
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